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f^) Abstract. We outline a remarkably efficient method for generating solutions to 

' quantum anharmonic oscillators with an x 2M potential. We solve the Schroedinger 

bJ) 1 equation in terms of a free parameter which is then tuned to give the correct boundary 

condition by generating a power series expansion of the wavefunction in x and applying 
a modified Borel resummation technique to obtain the large x behaviour. The process 
Q\ , allows us to calculate energy eigenvalues to an arbitrary level of accuracy. High degrees 

of precision are achieved even with modest computing power. Our technique extends 
to all levels of excitation and produces the correct solution to the double well oscillators 
\ even though they are dominated by non-perturbative effects. 
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t^- ■ 1. Introduction 

o : 

^ Harmonic oscillators are a corner-stone of many branches of physics. Consequently a 

large variety of methods have been used to study the eigenvalue properties of anharmonic 
oscillators (see pQ [2] and references therein for a general review) . High levels of accuracy 
have always been difficult to achieve due to slow convergence or often non- convergence 
of asymptotic perturbative expansions. For example the Bender Wu [3] expansion of 
| the quartic anharmonic oscillator ground state energy eigenvalue in positive powers 
of the coupling is known to be divergent for all non-zero values of the coupling. 
Methods of resumming asymptotic series [4] have been applied to generate approximate 
eigenvalues [5] [6j. In addition some types of anharmonic oscillators are dominated by 
non-perturbative effects such as instantons [7j. More innovative approaches have been 
required to produce a greater level of accuracy and account for these non-perturbative 
effects [8j[9|[10j[llJ. In addition to the numerical approaches some progress has been 
made in determining the analytic structure of certain anharmonic oscillators [12] [13] • In 
particular [13] outlines a type of anharmonic oscillator which is quasi exactly solvable 
with certain parts of the spectrum known exactly. 

These problems extend into quantum field theory. For example, the renormalisation 
group implies that the energy eigenvalues in Yang-Mills theory cannot be solved for 
perturbatively. Strongly coupled theories in particular are hard to deal with using 
traditional techniques. Anharmonic oscillators are therefore of great interest because 
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of their applicability in many branches of physics and because their mathematical 
properties often mirror those of other physical systems. 

We will outline a method of constructing solutions to the Schroedinger equation for 
an anharmonic oscillator of the form 

_ + px ^ + gx ^m = E m (i) 

dx 2 

lim ^ = (2) 

where x is real and units are defined to absorb Plank's constant and the mass such that 
h = 2m = 1. We do this initially by constructing a solution to the differential equation 
(CQ) in terms of one free parameter for a given p and g. We then vary this parameter 
until we observe the correct large x behaviour determined by the boundary condition 
(12]) using a contour integral method of resummation. We find the energy eigenvalues 
with an arbitrary level of accuracy. The process is easily automated to produce very 
high levels of precision even with modest computing power. 

In section [2] we will outline the basic method for the ground state of the quartic 
oscillator, M — 2. In section [3] we will extend the method to produce excited 
wavefunctions and energy eigenvalues. Finally in sectionH]we will show how this method 
can be extended to general anharmonic oscillators with an x 2M potential as in (CQ). 



2. Tuning for Large x 

In this section we will find the ground state wavefunction and energy eigenvalues 
corresponding to the quartic anharmonic oscillator obtained from (CQ) by setting M = 2. 
Since the ground state has no nodes we will construct the wavefunction in the form 
^ = exp(iy). We will make an even powered x expansion W = Y^Li a nX 2n since 
both the potential term and boundary condition are even. The coefficients a n are then 
determined in terms of the parameters p, g and E via (JTJ). Having chosen two of these 
parameters the third must be determined by ensuring the correct boundary condition 
([2]), which implies that W ~ — v /^x 3 /3 for large positive real x. Since our expansion 
for W in positive powers of x is only valid for small x we shall resum by analytically 
continuing x into the complex s = 1/x plane and using Cauchy's theorem to examine 
the large x behaviour. We define 

where C is a large circular contour about the origin in the complex s plane. The large x 
asymptotic behaviour implied by the differential equation requires W(s) to have a third 
order pole at the origin. This contributes a term —y/g/18 to L(X) by Cauchy's theorem. 
When the boundary condition is satisfied we find that any remaining singularities of 
W(s) lie to the left of the imaginary axis. The contribution from these is exponentially 
suppressed in L(X) so that in the large A limit only the singular contribution at the 
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origin remains, lirn,\-»oo L(X) = —y/g/18. In reality ([3]) is not calculated exactly but by 
truncating W at some order x 2N . Thus 

1 1 r f Xs N n N \2n-3 

JC n=l n=l v ' 

where in the evaluation of the contour integral we used the identity J c ds s~ n exp(As) = 
27ii\ n - 1 /T(n) for n < [H]. 

We proceed by finding our x expansion in W and look for the correct behaviour 
in Lat(A). We will consider solutions with a fixed coupling g = 1 and look at the 
relationship between E and p. We do this without loss of generality since our parameters 
are related by scaling properties of the Hamiltonian, as first noted by Symanzik and 
discussed in [15]. To help us we will scale x — > cx (c G R) in the differential equation 
(Cp) in such a way that we are free to place a restriction on our expansion for W . We 
can choose k = ai/a 2 at least up to a sign, say k = ±4. Now substituting W into our 
scaled differential equation and comparing coefficients of x 2n 

Ec 2 = -2ai (5) 
pc 4 = 4a? + 12a 2 (6) 
c 6 = 16aia 2 + 30a 3 . (7) 

We eliminate c to find expressions for E and p in terms of a 2 and 03 

^ —2kao , \ 

E = r 8 

(16^ + 3003)3 

Ak 2 a 2 2 + 12a 2 

(16A;ai + 30a3)i 

whilst for n > 3 we have 

a n+1 = - I ^ 4m(n - m + l)a m a n _ m+ i J / (2(n + l)(2n + 1)) , (10) 

\m=l / 

giving a n+ i in terms of a 2 and 03. 

Our goal is now to determine 03 for a given a 2 in such a way that the boundary 
condition is satisfied. We do this by tuning a 3 until the correct large A behaviour is 
observed in Ltv(A). To illustrate the process we shall choose positive k, k = 4 and 
a 2 = —3/16. With this sign choice and a 2 we get a zero p term. We choose a fairly 
modest iV initially, guess a value of a 3 then plot L N (X) and L N _i(X). L N and L N _i 
only provide a good approximation to £(A) for values of A up to the point where they 
appreciably diverge from each other. Therefore we restrict our consideration of A to 
within this range. 

With 03 too small we encounter a curve rapidly decreasing such as in figure [Th. 
With a 3 too large we encounter a curve rapidly increasing as in figure [Tb. An optimal 
value of 03 will give a curve flattening as we increase A as in [lb. We tune 03 until we 
achieve this. As a 3 gets closer to its correct value the exponential behaviour in figures 
and db becomes less pronounced within our range of acceptable A and flatness becomes 



Solving the Anharmonic Oscillator: Tuning the Boundary Condition 4 




a) a 3 = 0.015 - too small b) a 3 = 0.0193604 - optimal c) a 3 = 0.025 - too big 
Figure 1. L N (X) with N = 19,20 for a 2 = -3/16 




a) a 3 as a function of 02 b) p as a function of a 2 c) E as a function of 02 



Figure 2. p, J5 and 03 as functions of 0,2- 

a less well defined concept. This determines our level of accuracy for determining 03. To 
improve our accuracy we must increase N in order to consider larger A. As we consider 
these larger A we again encounter the exponentially increasing / decreasing behaviour 
which enables us to further tune a 3 to a greater accuracy. 

We completed this procedure in this zero p case and determined 03 to 6 significant 
figures with N = 20. With N = 100 we tune a% to 30 significant figures and with 
N = 300 we get 

a 3 = 0.01936043720245950419201997531721233596425589581549397570027615152, (11) 

to 65 significant figures and we find L(\) ~ —0.0934774 which is remarkably close to 
the predicted value of —0.0934723. This calculated figure of a 3 is accurate to the stated 
number of digits i.e. 65 significant figures and in agreement with existing literature 
[8j[9j[10j[llJ at least up to the 10-16 significant figures they quote. With a 3 determined 
we calculate the ground state energy via (jSJ) 

E = 1.0603620904841828996470460166926635455152087285289779332162452417 (12) 

again quoted accurately up to 65 significant figures. 

Calculating large numbers of terms is easy, even with modest computing power, 
given the linear nature of the calculations. The tuning process is easily automated. 

We repeat this for various a 2 and plot the results in figures [5] and El The two 
branches correspond to differing sign choices of k. With k = +4 we found solutions 
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Figure 3. £ as a function of p 



corresponding to positive energy. The solutions have a positive p term for a 2 < —3/16 
and a negative term for —3/16 < a 2 < 0. With k = —4 and > a 2 > —3/16 we found 
negative energy eigenvalues corresponding to p < 0. 

We also verify that non-zero p terms correspond to the literature by for example 
calculating the p = — 1, g = 1 energy eigenvalue. In doing so we must tune a 2 with 
k = +4 to obtain the correct p term. We found that a 2 = 0.004048768355681543705 
approximated p = — 1 with an error in the order of 10~ 16 . This value of a 2 is within ±5~ 21 
of the correct a 2 required to evaluate p exactly. The energy eigenvalue produced from 
this approximate value of a 2 gave us the same eigenvalue as stated in previous literature 
to within the 16 significant figures available for comparison. This is an example of an 
eigenvalue where instanton effects would normally dominate and perturbative techniques 
in H or g would fail. 

We now explain why this method of tuning is so sensitive. With M = 2 the 
differential equation (JTJ without the boundary condition (j2J) in general has an asymptotic 
large positive x solution of the form 



For A < 0, i&i has zeros along the real x axis however for A > 0, has zeros in 
the complex x plane off the real axis. Our boundary condition (J2J) requires us to take 
A = in which case has no zeros. We note that for A / 0, log will have a 
pole (possibly part of a cut). Such a pole contribution in the right half x plane would 
spoil our resummation of the large x behaviour. We have numerically determined the 
location of zeros in our wavefunctions for varying a 3 and shown that they numerically 
approximate the location of the zeros in our asymptotic large x solution for varying 
A. Thus varying a 3 corresponds to varying A in ( 113ft . The presence of these poles is 
responsible for the rapidly increasing /decreasing behaviour for values of 03 on either 




(13) 



Solving the Anharmonic Oscillator: Tuning the Boundary Condition 



6 



side of the correct one due to the exponential factor in ([3]). It is this behaviour that 
allows us to select the correct value of 03 to any specified level of accuracy. 

3. Excited States 

We now construct the excited states and energy eigenvalues of the quartic anharmonic 
oscillator. Firstly we write the qth excited state as \& q = P q ^o where the energy is 
E = E q + E , and \l/ is the ground state obtained in the previous section. For q odd 
P q is odd and for q even P q is even. We therefore expand P = J2^=o c ™ x ™ anc ^ sum on ^ 
over even or odd values of n as appropriate. We set either c or c\ to unity as a choice of 
normalisation. The remaining c n and E q are then solved for using a recurrence relation 
in terms of either c 2 or C3. This is easily obtained from our new differential equation 
which comes from substituting ^> q into (CQ) to obtain 

d 2 P + 2 dW dP q + E p _ Q 

doc doc doc 

This differential equation has two types of large x solution. Either 

P~exp( or P ~ exp (2^/g~x 3 ) . (15) 

V 2y/gxJ 

For the correct boundary condition (j2J) we must choose the first type of solution. We 
therefore construct 

1 r f> Xs N r - \ 2na 

TjV ( A ) = _L / ds — Y^ = Yc n —^ (16) 

v ; 2m Jr s s 2na ^ T(2na + 1) v ; 

Ji -' n=0 n=0 v ' 

and look for a flat curve as we tune c<i or C3. We have introduced an additional 
parameter a by substituting s — > s a in P(s) since we find that P(s) has a more limited 
region of analyticity than W(s) when the boundary condition is satisfied. Here we only 
assume that P(s) is analytic in some wedge shaped region radiating from the origin and 
containing the real axis. Singularities outside of this region of analyticity are observed 
in T/v(A) in the form of oscillations. They can however be rotated in the complex s 
plane so that they lie to the left of the imaginary axis by reducing the parameter a < 1 . 
Having done this the singularities become exponentially suppressed. 

We illustrate the process in the zero p case for the odd eigenf unctions. There will be 
multiple values of r = — C3 that correspond to different levels of odd excitation. Let us 
label these r n in such a way that r n+1 > r n . With r < T\ we obtain a rapidly increasing 
curve however with T\ < r < r 2 we get a rapidly decreasing curve (figure HJ) . We follow 
our tuning procedure in the same manner as for the ground state however this time we 
do not encounter a flat curve but oscillations. These result from a pole or cut outside of 
our region of analyticity. We could however take a smaller a to recover a flat curve and 
proceed with our tuning procedure. For r = T\ we found a = 0.6 sufficient to achieve 
this. 

We can produce the full spectrum of eigenvalues by continuing to vary r. We 
find that as r passes through a value r n we switch from the rapidly growing to 
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d) t « t 2 « 1.99546 e) r = 2 

Figure 4. T(A) with an odd prefactor 



rapidly decreasing behaviour. With r 3 > r > r 2 for example we switch back to the 
rapidly increasing curve. This alternating behaviour continues with higher excitations 
as illustrated in figure HI Exactly the same procedure works for even excitations but 
we vary c 2 instead of C3. Having found an eigenstate through this method we cannot 
immediately tell which energy level it corresponds to. To do this we could plot the 
prefactor using a similar contour integral method of resummation. We then count the 
number of nodes. We did this for some of the lower excitations. We calculated excited 
states up to q = 39 with g — 1 and again found exact agreement to the quoted level of 
accuracy in previous literature [8][9j[10j[llj. We give some of these eigenvalues in the 
appendix. 

Whilst we cannot attribute the rapidly increasing / decreasing behaviour of T/v(A) 
to zeros in P q we believe that a similar effect is encountered this time due to the large x 
behaviour. There were two types of large x behaviour (|15|) that we were able to derive 
from the differential equation (JHJ). We chose the first in order to satisfy our boundary 
condition When c 2 or c 3 do not correspond to an energy eigenstate we believe that 
we are obtaining the second type of solution. We note that such large x behaviour would 
give an additional pole contribution to T(A). Again our resummation is conveniently 
spoilt. We have numerically verified this result by plotting P for large real values of x 
for a range of C3 by exploiting Cauchy's theorem. 
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4. Other Potentials 

In this section we consider other values of M in [U The large positive x behaviour is 
now W ~ —^gx M+1 /(M + 1). We should therefore redefine our L^v(A) for a general 
x 2M potential 

1 1 f e Xs - a n N \(2n-M-l)a 

L " {X) = 2^W^ c ds — Y.^ = Y. a - T{2na + iy ^ 

JC n=l n=l v 7 

where again we introduce the parameter a since for M > 2 we find that W(s) is 
analytic within a more limited region. Our prescription of reducing a < 1 will therefore 
be required to rotate these singularities to the left of the imaginary axis where they 
become exponentially suppressed. 

Our a n are again determined via the differential equation in the same manner as 
before. We apply the rescaling x — > cx so that we can fix ai/a 2 = ±4 as before. We 
pick a value of a 2 and use (fTUl) to solve for all of the coefficients in terms of a n+ i. This 
relation now holds for n > 2 but not n = M. In its place we have 

M 

c 6 g = 2(M + 1)(2M + l)a M +i + ^4n(M -n+ l)a M a n -M+i (18) 

n=l 

which is then substituted into (j3J) and (jSJ) to give E and p. 

Our procedure is now the same as for the M = 2 case. We do find however that 
for increasing M the region of analyticity becomes smaller and therefore an increasingly 
small a is required. We performed this procedure with M ranging from 2 to 50 for g = 1 
and found results matching those in [9] for M = 2,3,4. Having determined the ground 
state we have applied the technique outlined in section [3] to obtain some excited energy 
eigenvalues. Again these are in complete agreement with [9]. 

5. Summary 

We have developed a method for calculating the relationship between the physical 
parameters of a general x 2M anharmonic oscillator. The equations we solve are linear 
and the process of refining our estimate is easily automated. We can calculate the 
physical quantities and wavefunctions for all levels of excitation to an arbitrary level of 
accuracy with an error that can be reduced by increasing the number of terms in our 
expansion. Using modest computing power we have demonstrated that high degrees of 
accuracy can be obtained very quickly. Our technique overcomes some of the deficiencies 
of traditional perturbative techniques which rely on coupling constant expansions and so 
do not immediately reveal the effects of instantons, for example. Finally, we note that 
the analytic continuation of quantum mechanical systems into complex configuration 
space has recently been studied in VT symmetric quantum mechanics (see [TU] and 
references therein). We believe that understanding the properties of Hermitian theories 
in the complex plane is still of great interest. 
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Finally we note that in the case of the quasi exactly solvable solutions studied in 
|13j . the expansion of both W and P in powers of x becomes truncated. In this type 
of solution it is more obvious that the correct boundary condition is satisfied by the 
large x behaviour. This is trivially reflected in our resummation technique. We have 
numerically verified that the results of p2] are correctly reproduced for some specifc 
choices of an x 6 polynomial potential. 

Appendix 

Below we give some of the excited energy eigenvalues for the p = quartic (M = 2) 
anharmonic oscillator. The results represent accurate eigenvalues rounded to 48 
significant figures. 

q E Q + E 1 

1 3.79967302980139416878309418851256895776606546733 

2 7.455697937986738392156591347185767488137819536750 

3 11.6447455113781620208503732813709364365508721620 

4 16.2618260188502259378949544303846135342445865045 

5 21 .23837291823594002414971 1 1 135886363767048320597 
20 122.604639000999455020762971417615181874976633223 

38 284.068590581400743150496281208125064777084713267 

39 293.948458266006085433669997483521626303445899275 
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